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Nomenclature 


A area 

Bi Biot number, eqn (2-6) 

c specific heat 

D domain 

Fo Fourier modulus, eqn (2-5) 

h convection coefficient 

Bessel function of the first kind of order t 

k thermal conductivity 

L length of cylinder 

m mass 

P point 

Q heat transfer 

^ total heat transfer 
tot 

q heat transfer rate 

qt, ■ constant heat flux 
R radius of cylinder 

r radial coordinate 

T temperature 

T average temperature 

t time 

V ' volume 
z axial coordinate 

a thermal diffusivity _ 

3 < 


"wr 





6 Kronecker delta 

cp angular coordinate 

$ size parameter 

p density 

Subscripts 
c center 

i initial or evaluated at node 

j e valued at node j 

s surface 

w surroundings 



1.0 INTRODUCTION 

In the processing or coohing of foods, nutrient materials are 
heated up to the sterilization temperature and "retained at that 
temperature for sufficient time to destroy harmful bacteria. When 
foods are cooked immediately prior to serving, the nutrient material 
simply cools from temperatures in the sterilization region down to 
serving temperatures 135-150F. When foods have been processed pre- 
viously and need only to be heated for serving, the nutrient material 
is heated from some initial temperature (whether room temperature or 
refrigerated temperature) up to serving temperatures. The cool, pro- 
cessed foods have a low level content of harmful bacteria. However, 
the rate of bacteria growth is accelerated in the temperature range 
45-140F. Therefore, it is important in the heating process for 
serving to be accomplished in a reasonable time so that nutrient 
material does not remain in the region of accelerated bacteria growth 
for too long a time. For example, it is possible for bacteria to 
double by cell reproduction in just 20-30 minutes. Therefore, in 
the heating process foods cannot linger in the critical temperature 
range. 

For the heating of foods in the ordinary. Earth-based facility, 
the primary mode of heat transfer is the convection mechanism, which 
is a very effective mechanism of heat transfer for substances having 
fluid characteristics. Since water boils at 212F, the surfaces of 
the container can be maintained at temperatures well above that of 
the food substances without permitting boiling of the food. The 
presence of large temperature differences increases the rate of heat 
transfer. Convection heat transfer and the presence of large temper- 
ature differences each enhance the heat transfer and reduce the time 
required for the heating process in the Earth-based facility. 

In the space vehicle, the heating of foods is desirable for the 
comfort and for the psychological benefits of the crew. However, 
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the heat transfer is impaired by two physical factors; (a) the 
zero-q envir^^#lfent and (b) the reduced pressure level of approxi- 
;.mate'ly 5 psia. Since the convective mode of heat transfer depends 
upon bouyant forces, the zero-g environment eliminates convection; 
the food material must distribute energy internally by the conduc- 
tion mode, which is a less effective mechanism of heat transfer . 
Also, at the reduced pressure level, water boils at 160F . Therefore 
the temperature of the walls of the container cannot be elevated 
substantially above the temperature of the food material — thus, 
eliminating the beneficial large temperature difference. 

There are then two primary tasks to be accomplished through 
the thermal modeling of nutrient systems: (a) prediction of the 

time required to heat foods to desired temperatures, and (b) devel- 
opment of parametric studies to optimize the system for minimum 
power consumption. 

1-1 Food Heating System for SkyLab 

The physical system for the heating of foods for SkyLab is a 
tray arrangement with several receptical cavities for the insertion 
of canned foods. The system uses two sizes of aluminum cans. The 
cavities are lined internally with blanket-type, electrical resis- 
tance heaters. The heaters are thermally controlled and provide a 
uniform heat flux of 2 watts per square inch. The sensor for the 
heater control system is a thermocouple, which is attached to the 
cavity wall. 

The control system turns the heater off when the temperature 
sensor reaches 155F and reactivates the heater when the temperature 
drops to 143F. The food heating system is designed to provide hot 
foods at 149±6F. The aluminum cans of food receive heat on the 
sides and on the bottom with the top insulated. There are several 
different initial states of the food when inserted into the heating 



. system. The contents of the cans may be initially at a uniform 

« 

temperature of -lOF (frozen storage), 60F (ambient storage) or 
-.ri30F'"(rehydrated) . The contents of the can are pressurized to 
5 psia with nitrogen. The requirements of the food heating system 
is that foods be heated above 140F within 132 minutes; this time 
requirement assures that the food passes through the microorganism 
growth temperature range (45-140F) rapidly enough to prevent con- 


tamination . 
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. 270 BASIC THERMAL MODELING OF NUTRIENT MATERIALS 

i 

-The theiilfuaC analysis of the heating or cooling of nutrient 
“materials can be considered in two parts. First, the establish- 
ment of the mechanism (or mechanisms) by which heat is transferred 
to (or from) the nutrient material; and second, the determination 
of how heat is transferred within the material itself. 

The basic mechanism by which heat is transferred (or from ) 
the material can involve any combination of the three basic heat 
transfer mechanisms (conduction, convection, radiation) . Heat 
transfer to the material by conduction involves the direct contact 
of two solids where contact resistance may be a consideration. 
Convection describes the heat transfer mechanism between a solid 
surface and a fluid. While thermal radiation is always present, 
it is particularly important in the absence of a transferring media 
since the other modes are non-existarit . However, in order for ther- 
mal radiation to be significant, relatively large temperature dif- 
ferences must be present. 

Two basic models are commonly involved in the analysis of the 
transfer of heat within the nutrient material, itself. The nutrients 
may be solid or plastic so that temperature gradients may be present 
in the material. On the other hand, the nutrient material may have 
sufficient fluid character so that moderate mixing occurs (whether 
by movement within the container or by natural convection currents). 

In this latter case, the temperature throughout the material is 
almost uniform (for moderate heating rates) .• 

2.1 Nutrient Materials with Mixing 

When a nutrient material has some fluid characteristics so that 
mixing may occur, the presence of internal temperature gradients are 
minimized, and the system can be characterized as having a uniform 
temperature at any instant. This basic model (illustrated in Fig. 2.1) 
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is known as Newtonian heating (or cooling) and is the earliest 

( 

model employed in transient heat transfer analysis. At any par- 
:^ticTul'ar instant, the nutrient system has a uniform temperature. 



S 


Fig. 2.1: Newtonian Heating Model 

An energy balance on the system provides 

dT „dT 

q = mcr— = pcV— 


(2-1) 


If the heat transfer at the boundary is by convection, the heat 
transfer rate can be expressed in terms of the convective mechanism 
q = hA [T„ - T(t)] (2-2) 

If eqns (2-1) and (2-2) are combined and integrated, the tempera- 
ture response of the system becomes 


Too - T(t) _ k V 

- T, ~ ® . 

1 

For a cylinder of radius R and length L 


A 

s 

V 


2/R 
— ( — + 

R\L 



(2-3) 


(2-4) 


If the Fourier modulus is defined as 



(2-5) 
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’and the Biot number is defined as 

< 



then, the dimensionless temperature response of the system can be 
written as 

= exp[-2 BiFo (l + 0] (2-7) 

which is shown in Fig. 2.2. 

The Newtonian model also applies to solid materials whose 
thermal conductivity is very large (such as silver or copper) ? 
a more accurate statement would be that the internal resistance 
to heat transfer by conduction is very small as compared to the 
surface resistance to heat transfer by convection (Bi < 0.1). 

The Newtonian model relations can be modified to include the thermal 
radiation mechanism; however, the ensuing integration is not straight- 
forward. 

2 . 2 Nutrient Materials Without Mixing for Fixed-Temperature Boundary 
Conditions 


When a nutrient material is of rigid texture (i.e., sufficiently 
solid to avoid mixing or convection currents) , the transfer, of energy 
internally is by the conduction mechanism. The basic partial differ- 
ential equation for the axisymmetric cylindrical configuration (as 
shown in Fig. 2.3) may be written as 


^r^ r ^r a St 


(2-8) 


for the case of constant thermo physical properties. The solution 
of eqn (2-8), subject to specified boundary and initial conditions, 
provides the local temperature distribution T(r,z,t) at any instant. 
From the temperature distribution, the response of the average tem- 
perature of the material (or the center temperature) to a sudden 
change in environment at the surfaces of the finite cylinder can be 
established. 






Fig. 2.3: Analytical Model of a Finite Cylinder 

Analytical solutions of eqn (2~8) are available [l,2]* for 

basic boundary and initial conditions. For example, if the nutrient 

material is initially at a uniform constant temperature T. 

1 

T(r,z,0) = (2-9) 

and if the surfaces of the cylinder are suddenly changed and main- 
tained at a uniform temperature 

T(r,L/2,t) = T(r,-L/2,t) - T(R,z,t) = T (2-10) 

s 

The solution of eqn (2-8) for the initial and boundary conditions 
(2-9) and (2-10) is 

♦Numbers in brackets indicate references. 
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TClTfZjt)” T t 0*7\ ^ \ T 

-T - T Z Z R (^m ^ n U^ J 

i s n=l m= 1 *• 


( 2 - 11 ) 


where 


A = 


n (Un) 


and 


^ (“1) 
m 


m+1 2 




n 


S (2m - 1)~ and Jq (wtn) = 0 


The evaluation of eqn (2-11) involves the summation of a Fourier 
series which converges very slowly; special techniques are generally 
required to accelerate the convergence. From the evaluation of 
eqn (2-11) at a particular point, the response of the temperature 
of that point with time can be determined. For example, if the 
surfaces of a cool cylindrical container at a uniform temperature 
were suddenly changed to a warmer temperature T^, it would be of 
interest to know how fast the "cold spot" in the system responded. 
Equation (2-11) would then be evaluated at the centroid of the 
homogeneous cylinder (r=0, z=0) . Fig. -2.4 shows the response of 
the dimensionless center temperature, T^, with dimensionless time. 

It is also of interest to know how the average temperature of 
the system is responding which can be obtained by 


1 r'" 

T(t) - “ I T(r,z,t)dV 
V Jo 


( 2 - 12 ) 


where eqn (2-11) provides the temperature expression for substitu- 
tion into eqn (2-12) . 


T(t) - 

T. - T 
JL S 


■„li 1.1 k k ■ 


(2-13) 


Fig. 2.5 shows the response of the average temperature of a cylinder 
in dimensionless form. 

Additional information can be gained from the solution of 
eqn (2-11) . The amount of heat transfer can be established by any 
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Fig. 2.4: Response of the 
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of several methods. Of particular interest is the fraction of the 
total possible heat transfer which has occurred up to any particular 
time. . For example, if a long enough time is provided the entire sys- 
tem will assume the surface temperature T and the total heat . transfer 

s 

to the cylinder can easily be determined by the change in the internal 
energy of the cylinder, namely 

Q = mc(T - T.) (2-14) 

tot S X 

The heat transfer to the cylinder from t = 0 to t = t can be 
written as 

. Q(t) = mc[T(t) - T^] (2-15) 

The fraction of the total heat transfer transferred to the cylinder 
during the time interval t = 0 to t = t can be expressed in terms 
of the average temperature of the cylinder 


Q(t) 

^tot 


T,. V - T. 
(t) 1 

T - T. 
s X 


(2-16) 


Since the right-hand side of eqn (2-16) is presented graphically 
in Fig. 2.5, the dimensionless heat transfer- is also represented 
by Fig. 2.5. The fraction of the total heat transfer is also pre- 
sented in different form in Fig. 2.6. 

"Fig. 2.2 through 2.6 present the thermal characteristics of 
finite cylinders for various size cylinders (i.e., for various L/R 
ratios) . For the cases presented, the cylinder is exposed to iden- 
tical conditions on each surface. (i.e., the radial surface and the 
two ends of the cylinder). However, the results can be used for 
the finite cylinder, which is insulated on the ends; this case cor- 
responds to an infinite L/R ratio where end effects are neglected. 

2.3 Modeling Consideration for SkyLab Configuration 

As discussed in the Introduction, the absence of bouyant forces 
eliminates natural convection currents in fluid-lihe substances. The 
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L/R 







distribution of heat within the nutrient material is then dependent 
ion the conduction mode. The conduction mechanism can be adequately 
described by energy relationships. The two-dimensional, transient 
conduction equation in cylindrical coordinates was given as eqn (2-8) 
for the case of constant, uniform thermophysical properties. The 
solution of the partial' differential equation for a particular appli- 
cation requires the establishment of appropriate initial and boundary 
conditions which simulate the physical situation. 

There are unlimited combinations of possible boundary and ini- 
tial conditions. However, the many possibilities generally can be 
categorized into one of three broad classifications: (a) specified 

surface temperature; (b) specified surface conductance; and (c) 
specified surface heat flux. In the proceeding section, the surface 
temperatures were specified; solutio-ns associated with this type of 
boundary condition are usually not too difficult, but, the boundary 
condition does not accurately describe real-world situations. The 
other two classifications of boundary conditions occur more frequently 
in nature. The specified surface conductance describes the inter- 
action of a fluid and a solid surface through the convection mechanism. 
Although the surface conductance (or convection coefficient) is 
assumed constant, the heat transfer rate is not if there is a change 
in the temperature of the solid surface and/or the fluid. The third 
classification of boundary conditions is the specified surface heat 
flux, which, accurately describe situations involving electrical 
heating elements, such as in the ShyLab system. In the SkyLab con- 
figuration, a uniform heat flux acts on the sides and bottom of the 
container with the top insulated. In order to simulate the SkyLab 
system, the model must be capable of describing intermittent heating 
periods. 

Due to the presence of zero-g conditions, the boundary condi- 
tions may be complicated somewhat. If the nutrient material fails 
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. to "wet" the surface of the container, the food may "float" in a 

t 

partially filled container. The situation may occur where the 
_--food‘“separates from the heating surface. The nitrogen gas filling 
the voids surrounding the nutrient material has a very low thermal 
conductivity and acts as an insulator between the heating surface 
and the nutrient material. In this situation, thermal radiation 
may be the predominate mechanism to be considered. Since the tem- 
perature of the heating surface is only 149±6F, the radiation mech- 
anism is not an effective mode of. heat transfer. Since heat is not 
effectively transferred to the food, the heater is only heating the 
wall of the can. When its temperature exceeds the cutoff, the heater 
is deactivated. The heater will in fact be off most of the time, and 
the heating time will be extended substantially. 

In modeling food substances, the nature of the nutrient materials 
themselves is an important consideration. Foods have hetrogeneous 
character (i.e., a non-uniform composition such as soups, stew, etc.). 
It is standard practice in modeling such systems to establish some 
"weighted" average values of the thermbphysical properties for the 
various food components and then to treat the material as homogeneous 
in nature. In the zero-g environment, the nutrient material may 
form internal voids or cavities which may be filled with gas. The 
presence of voids in the material may effect the internal heat 
transfer substantially (as compared to the hetrogeneous character- 
istics of several food substances) because of the insulating effect 
of such cavities. The available data concerning the thermophysical 
properties of nutrient materials are somewhat limited. The data 
are scattered in many diverse sources in the literature- Also, in 
gathering data from the literature, one questions the validity of 
some of the data because of the experimental techniques employed in 
the measurements. 


\ 
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3.0 ANALYTICAL APPROACH 

3^1 Boundary and Initial Conditions 


For the model illustrated in Fig. 2.3, the conduction relation 
for the heating of nutrient materials was given in eqn (2-8) for 
the case of constant thermophysical properties. To solve this rela- 
tion, appropriate initial and boundary conditions must be selected 
which best describe the physical configuration of the Shylab system. 
If the nutrient material is initially at a constant, uniform tempera- 
ture, the initial condition becomes 

T(r,z,0) = T^. 

The boundary condition for the insulated top is 


?,T(r,L/2,t) ^ 

az 


(3-2) 


For the intermittent heating of the sides and bottom, the appropriate 
boundary conditions can be written as 


^ (constant heat flux; 
;^T(R,z,t) _ >T(r,-L/2,t) ^ heater on) 




az 


(3-3) 

0 (insulated; heater off) 


Along with the physical boundary conditions is the stipulation that 
the temperature remain finite along the axis of the cylinder. 


3.2 General Solution 

There is no single analytical solution which will satisfy all 
these conditions simultaneously. However, a piece— wise solution 
can be formed by adding together a series of solutions, each valid 
for a short time increment (e.g., one heating period) . There are 
actually only two different solutions required. One for the time 
when the heater is on, and the other for when, the heater is off 
(which is really a special case of the first with qo = 0) . 



Olcer [31 developed a general solution (additional details in 
Appendix A), For a flux, qo , applied at the radial and bottom sur- 
faces, insulated top, and- an arbitrary initial temperature distribu- 
tion, T(r,z,0) = T^(r,z) (For the second and succeeding heating periods 
the initial temperature is non-uniform.), the transient axisymmetric 
temperature response is:- 

L . ■ ‘ 

T(r,z,t) = 


+ 


R^L 


IkL 




0 

0 - L 


(I - “T 


12/ 2k\R^ 2/jlo J 


16tt 

R®L 


Jo(Mmr) cos [— + ~ | exp (-ecX^^ t) 
2(1 + 6„„) Jo® (HR) 


I I 

m=0 n^ 


no 


L 
2 R 


[ J J T.(r,z)rdrdz 


_ L 0 
2 


^ J.(MmR)fo°}] 


kX |a 
mn'^m 


(3-4) 


where indicates the two cases of heater on and off 


X ® 

mn m \ L / 


(3-5) 


where u R s 0 is the mth root of 
m 


Jo' (M^jR) = 0 


(3-6) 


and the prime ( ^) denotes differentiation with respect to the argument 
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For the first heating period the initial temperature is 
uniform and 

T. (r, z) = T. (3-7) 


so that the first integral becomes equal to T^, and the other 
integration involving vanishes. 

For a given initial temperature distribution and surface heat 
flux, eqn (3-4) can be used to calculate the local temperature within 
the cylinder, T(r,z,t). The highest temperature will be at the bottom 
corner, point (R,0) . Therefore, the temperature T(R,0,t) is monitored 
until it reaches the predetermined maximum value. At this time, t^ , 
the entire temperature distribution, T(r,z,tj), is evaluated; this 
temperature distribution becomes T^(r,z) for the second phase of the 
first heating cycle. With q© = 0 (heater off) eqn (3-4) is used to 
determine the temperature within the food. When the monitored temper- 
ature at the bottom corner (R,0) drops to the predetermined minimum 
value, the heater is again activated. The temperature distribution 
at this time becomes the initial condition for the first phase of the 
second heating cycle. Equation (3-4) with non-zero boundary heat flux 
is used to determine the temperature. This procedure is continued 
until the food is heated to the desired level. The average food tem- 
perature can^ be determined either by averaging the temperature distri- 
bution: 

L 

2 R 

T(t) = ’J J T(r, z, t)rdrdz (3-8) 

^kO 

2 


or by dividing the total heat added to the food by the heat capacity 
of the food. Hence 


(t) = qo(2TTRL + TTR^) tpri 


where t is the total "on-time" for the heater, 
on 


(3-9) 
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suits for an Infinite Cylinder 

«i 

The ana.Xytik^.cal method outlined above has been programed and run 

— ■ ■- ■ ■ 

Zfor~the spv^^'.al case of an infinite cylinder (n = 0, L *+ “) heated 
on the cylPT^-dr ical surfaces only [4], An infinitely long cylinder is 
equivalent to a finite cylinder insulated on both ends. A discussion 
of the 'program and a simplified flow chart appear in Appendix B. 

Figure 3.1 depicts typical temperature distributions as a func- 
tion of non-dimensional radius, r/R (with heater at 142±8F) . Curve 3 
represents the distribution when the wall temperature first reaches 

T (i.e., the end of the first heating phase). Curves 1 and 2 are 
max 

the distributions at one-third and at two thirds of the total heating 
time of the first phase. The remaining odd numbered curves indicate 
the distributions at the end of each successive heating phase. The 
even numbered curves indicate the distributions at end of the succes- 
sive insulated phases (i.e,, just before the heater is reactivated) . 

3.4 Conclusions 

While satisfactory results were obtained with the analytical 
model, it was decided not to be continued to the finite cylinder case 
for the following reasons: 

1. A discontinuity always exists when switching from. one 
solution to the other due to the discontinuous manner 
in which the flux is switched on and off. The effect 
of this discontinuity is minor. 

2. Considerable computation time is necessary and the time 
required will increase considerably by the inclusion of 
the double summation required for the finite cylinder. 

3. A finite difference model appeared more attractive. 
Computation time was reduced and the finite difference 

• model is more versatile. The inclusion of different 
boundary conditions and a study of heterogeneous effects 
are reasonable extensions. 

23 < 
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4.0 FINITE DIFFERENCE APPROACH 

4 . 1 Introduction to Discrete Methods 

The ultimate goal of discrete methods is the reduction of a 
continuous system to an equivalent lumped-parameter system which 
is suitable for solution on a digital comput.or. The basic approx- 
imation involves the replacement of a continuous domain D by a matrix 
of. discrete points (nodes) within D, as shown in Fig, 4.1 for two 
dimensions. • 



Fig. 4.1. Discrete Approximation of a 

Continuous Two-Dimensional Domain 


Instead of developing a solution defined everywhere in D, only 
approximations are obtained at isolated points (nodes) , P. . , Inter 
mediate values, integrals, and derivatives may be obtained from this 
discrete solution by interpolation. This mathematical discretization 
replaces the derivatives by discrete approximations, (e.g., finite 
difference approximations) . 

At a given node, at a given time, the values of the coordinates 
define a unique location. For an axisymmetr ic cylindrical coordinate 
system, the set (r^, t^^^) define such a location where r^ is the 

25 < 
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■^r^Tal coordinate, is the axial coordinate, and t^^^ is the time 
after n time steps. Hence the dependent variable (temperature in* 
■this_-case) at this point is 


T = T (r . , z . , t 
1 J 


(n)% _ <ji(n) 

j 


(4-1) 


At a given time, the nodes (r . , z.), (r. z.), (r . , z. ) and 

ID 1-1 3 1 D“l 

^^i 1' 1^ define a subregion in the shape of an annulus of 

inner radius r. outer radius r,, and thickness z. - z. - as 
1-1 1 1 i-i 

depicted in Fig. 4.2. 



Fig. 4.2. Subregion Formed by Nodes 
(r.,z.)i (r. ,,z.)/ (r.,z. ,), and (r. ,,z. ,) 

i' j' * 1-1' j' * X J~1 1-1' j-1 


4.2 Explicit versus Implicit Models 

The standard finite difference approximations for partial derive 

tives can be found in most applied mathematics texts [5,6,7] and are 

26 < 



•dealt with in more detail in the next section. In general, there 

t 

are two_ techniques available for . approximating the differential 
_-equa~t"ibn in the transient problem. These are the implicit and 
explicit methods. 

In the explicit technique, the future temperature at some node 

(i,j) i's expressed in terras of its present temperature and that of 

surrounding nodes; that is, (for a first order approximation) 

m(ri+l) _ m{n) m(n) m(n) m(n)'\ /4_2^ 

’^i.j " ' ' 


where the form of the functional relationship depends on the govern- 
ing differential equation. Since the initial temperature distribu- 
tion is known, this technique can be used successively at each node 
to advance the solution one time step (i.e., eqn (4-2) is applied k 
times if there are k nodes) . The process can be repeated indefinitely 
until the required time has elapsed. The accuracy of the solution 
depends largely on the spatial nodal spacing. The closer the nodes, 
(hence, a larger number of nodes for a given physical system) the 
more nearly accurate the approximation'. The fundamental limitation 
of the explicit representation is that for a given nodal spacing the 
maximum rate at which the solution can be propagated (in time) is set. 
(If this rate is exceeded the solution becomes unstable) . In general, 
the more sparse the nodes, the smaller the maximum allowable time step 
(t^n+1) _ At) becomes. This can be a severe limitation to the 

explicit method' for a given application. To increase the maximum At 
requires more nodes, which in turn requires additional computation 
at each time step. Thus, even though fewer time steps may be needed, 
each one becomes more tedious. General rules have been developed for 
the maximum time step allowable under various conditions [8]. 

The implicit technique removes the possibility of this instability 
at the expense of increased computational . complexity . In the implicit 
technique, the future temperatures at a given node (i,j) and the 
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surrounding nodes are expressed in terms of the present temperature 
of node (i,j); that is, (for a first order approximation) 


f (T 


( n+1) 

i, j-1' 


n+1 ) ^ 
i# j+1' 


rp ( n+1 ) 

i-l.j 


T(n+1) = ^(n) 

i+l,j' i,j i,j 


(4-3) 


where again the functional form of f depends on the specific problem. 
This relationship must be determined at each node leading to a set 
of simultaneous, linear algebraic equations to be solved at each 
time step. This method is unconditionally stable in time but, of 
course, can introduce large round-off errors if the time step is 
excessive. Again it is noted that while the implicit method seems 
more attractive because of its stability characteristic, the require- 
ment in the implicit model of solving a set (one equation for each 
node) of simultaneous equation can hardly be compared to the solution 
of a set of independent equations required in the explicit model. 

4 . 3 Finite Difference Approximations 


The following are the finite difference approximations required 


^ 1 
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Central difference 
(4-4) approximation of 

the first derivative 

Central difference 
(4-5) approximation for the 
second derivative 


The forwaird difference 
(4-6) approximation for the 
first derivative 


The other required approximation follows. 

Due to the complexity of the implicit model, it was decided to 
use an explicit model for the problem first. The explicit finite 
difference approximation for eqn (2-8) around the node (i,j) is: 

2S< 
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(4-7) 


Equation (4-7) applies for any nodal spacing. However, it is 

convenient at this time to assign a uniform nodal network. At each 

z. there are N radial nodes. Node 1 is on the center line of the 
3 

cylinder; node N is on the circumf erencial surface. The nodes are 


equally spaced; i.e.. 


Ar = 


N - 1 * 


(4-8) 


There are M equally spaced' nodes in the z direction. Hence, 

L 


Az = 


M - 1 * 


(4-9) 


The time steps are also of equal length and are denoted At. With 
the nodal matrix thus defined eqn (4-7) becomes 
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(4-10) 


The required temperature is T. .^*^■*'1'. From eqn (4-10) 
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If 


g ( At ) 
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(A2)^ 
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= M 
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then 

and eqn (4-'ll) becomes 
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(4-12) 


Equation (4-12) is the general form of the algorithm required to 

generate the temperature at all interior points (i.e. , those not 

directly effected by the boundary) . 

At the center line, i = 1 and the coefficients of . 

111. 3 

become infinite in eqn (4-12), Howeyer, by symmetry 
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so that for 1 < j < M 
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At the outer radius for a constant heat flux boundary condition 


^'r=R = 


(4-15) 


The finite difference approximation to eqn (4-15) is (at i=N) 
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or 
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(4-16) 


Substitution of eqn (4-16) into eqn (4-12) yields for 1 < j < M 
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At the top surface 
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(4-18) 


The finite difference approximation to eqn (4-18) is: (at j=M) 
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Substitution of eqn (4-19) into eqn (4-12) yields for 1 < i < N 
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(4-20) 


At the bottom surface 
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The solution at the bottom surface becomes for 1 < i < N 
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At the top edge, i = N and j = M so that from eqn (4 
(4-19) 

= 2M „ + Tl - 2M - 2M 1 

N,M r N-1,M L ^ zj N,M 


ft.. fn^ 2N— 1 ,, 2 { Ar ) 

+ 2M + r* M cift 

z N,M-1 2{N-1) r k 


At the bottom edge i ^ N and j = 1 so that from eqn (4-17) 
(4-22) 
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At the bottom center '^q i '^2 1 eqn (4-23) 
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(4-27) 


4,4 Application of the Finite Difference Approximat ions 


A thirty-six node model (6 radial & 6 axial nodes) was selected 
for the homogeneous cylinder. Equations (4-12), (4-14), (4-17), (4-20), 
(4-23) , (4-24) , (4-25) , (4-26) and (4-27) were used in a Fortran IV 
G program. A simplified flow chart appears in Appendix C. 

A simplified analysis LS 3 was made to determine the maximum prop” 
agation speed of. the solution. From eqn (4-12) it is apparent that 
if the term 


1 - 2M - 2M 
r z 

is negative, then the larger , the smaller . This result 

r,j 

is physically unreasonable. Therefore, time increments must be re- 
stricted such that 

1 - 2M - 2M > 0 

r z • 
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For the can sizes expected oii Skylab (R < 2 inches, L < 4 inches), 

with a thermal diffusivity of water {'a ^ 5 x 10"^ ft^/hr) and for the 

R It ' 

nodal network assumed (Ar = -, Az = — ) , the restriction in eqn (4-28) 
becomes 


or 



1 

10“^ (At) 


25 25 100 

1/36 1/9 (At) 

100 

25(36)+ 9(25) = -089 hr ~ 5 minutes 


Hence the time steps cannot exceed 5 minutes each. From experi- 
mental evidence, it is known that as the can approaches its maximum 
temperature, the heating phase of heater cycle may be as short as a 
few seconds. If the model' is to be capable of simulating this heating 
period, the time step should be on the order of a second. Since the 
time steps required are well below the maximum allowable, the explicit 
model will be sufficient. 

The predictions of the finite-difference model can be demonstrated 
by looking at the graphical presentations in Fig. 4.3, 4.4 and 4.5. 

The temperature specification for the heated food is 149i6F. There- 
fore, the thermal control on the uniform heat flux heater are speci- 
fied accordingly. When the hottest point on the heating surfaces 
reaches 155F, the heater is deactivated and the food stuffs are treated 
as an adiabatic system (i.e., no heat losses); when the hottest spot 
on the surface has dropped to 143F, the heater is reactivated, and the 
cycle continues until the average temperature of the food reaches 149F. 
In each of the cases shown, the food was assumed to have thermal dif- 
fusivity of water and a uniform surface heat flux of 2.0 watts per 
square inch was employed. In each of the figures, two temperatures 
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are shown - the average temperature of the food and the temperature 
of the coldest internal point. 

The temperature responses of the small container and the large 
container are shown in Fig. 4.3 and 4.4, respectively. Two initial 
temperature conditions (T^ = 130F and T^ = 60F) are presented in 
each figure. This same model was employed in evaluating the tempera- 
ture responses for an initial temperature of -lOF and the results 
are shown for both container sizes in Fig. 4.5. It should be care- 
fully noted that this particular model did not have the capacity for 
including the latent energy associated with the phase change in the 
thawing process. If the nutrient materials were similar to water, 
the latent heat of fusion (or melting) is 79.7 calories per gram 
(143 BTU/lbm) . If the melting process were approximated by a "lumped” 
parameter (i.e., assuming that there are no temperature gradients 
present in the nutrient material during the melting process), the 
additional time required would be 

D.30 hours (small container) 

0.39 hours (large container) 

The finite-difference model can be modified to include the latent 
energy of thawing to provide more accurate predictions of the re- 
quired heating times. 

The basic finite-difference model has a number of input parameters. 
Each of these parameters can be investigated individually to determine 
its influence on the temperature response of the nutrient material. 


35 < 
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5.0 BASIC PARAMETRIC STUDIES 

The heat transfer analysis provides the vehicle for perform- 
ing parametric investigations. Parametric studies involve the 
variation of each physical parameter individually to establish 
the effect of that physical parameter on the thermal response of 
the system. Thermal considerations in the design of food ‘heating 
systems include the effects of the following physical quantities: 
(1) the thermophysical properties of the nutrient materials; (2) 
the power rating of the heater (i.e. , the energy per unit surface 
area) ; (3) the control temperatures which activate the heating 

element; (4) the initial temperature of the nutrient material; 
and (5) the dimensions of the container. In the heat transfer 
field, the temperature is generally presented as a temperature 
difference ratio and time in dimensionless forms using Fourier 
modulus; however, the following results employ real time and 
temperature. 

5 . 1 ; Effect of Thermophysical Properties 

The thermophysical properties of the nutrient material have 
a substantial influence on the requirec3 heating time. In the 
conduction mechanism, the food must be able to conduct energy from 
near the heating surface to the center of the container. Of 
particular interest is the thermal diffusivity q. which is the 
thermal conductivity divided by the product of the density and 
heat capacity (or specific heat) i.e., a = • Rather than 

showing the temperature response of particular nutrient materials, 
the results are shown in terms of a ratio of thermal dif f us ivities 
(more specifically, the thermal diffusivity is normalized by 
(dividing by the value for water at standard conditions) * The 
response of the average temperature of the nutrient material is 
presented in Fig. 5.1a for the smaller container and Fig. 5.1b 

‘S6< ’ 



Fig, 5.1a: Effect of Variation in Food Properties 

(Small Container) 
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for the larger container, 

' 5.2 Effect of Heater Output 

If a nutrient material is heated by a uniform heat source 
which is not controlled, the heater would remain activated 
continuously until the desired average temperature of the food 
had been attained. However, if the nutrient material near* the 
heating surfaces is not to experience excessive local temperatures, 
a control mechanism must be included in the heater circuit. The 
temperature-time relationships for the temperature controlled 
heater are compared to the uncontrolled heater in Fig, 5.2a for the 
small container and Fig. 5.2b for the larger container. The 
average and "coldspot" temperatures are depicted. 

If the heater were continuously activated, a decrease in the 
heater output by one-half would increase the heating time (to a 
particular temperature) by two. However, when a temperature- 
controlled heater is employed, this effect is altered substan- 
tially. The effect of a variation in the surface heat flux is 
demonstrated in Fig. 5.3 for the small container when the initial 
temperature is 70F. The interesting point is that when the uniform 
surface heat flux is reduced by a factor of 8 from two (2) 'watts 
per square inch down to one-fourth {h) watts per square inch, the 
heating time did not increase by a factor anywhere near 8. For 
example, the nutrient material was heated from 70F to 130F in 0.62 
hr for the high surface hfeat flux case and in 1.07 hr for the low 
surface heat flux case. In this situation when the surface heat 
flux was decreased by a factor of 8, the heating time was not even 
doubled. In optimizing a heating system for nutrient materials, 
the output of the surface heater is an important consideration. 

The reason the lower surface heat flux does not require a propor- 
tionally longer period of time is that the heater remains activated 
a substantially longer period of time prior to the initiation of a 

S9< 
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Fig, 5.2a: Comparison of Thermally Controlled Heater 

with a Constant Continuous Heater 
(Small Container) 
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periodic heater action. For the case of the higher surface heat 
flux, the heater begins its cycling pattern early. 

5.3 Effect of Temperature Controls 

The level of the temperature controls on the uniform surface 
heater has a great influence on the required heating times for the 
nutrient materials since the longer the heater • cycling process can 
be delayed, the shorter the required heating time. Although the 
Skylab system has established the temperature control range of the 
heater, an investigation of the control levels can provide insight 
into the operational characteristics of the system. 

In studying the effect of the thermal control of the surface 
heater, two temperatures are considered - the surface temperature 
at which the heater is cut off and the' surface temperature at which 
the heater is reactivated. In the model of the cylindrical con- 
tainer for the Skylab system, the point affected greatest by the 
surface heaters is the lower corner; the lower corner is the 
intersection of the bottom heated surface with the curved-side 
heated surface which corresponds to the location r = R, z = 0 in 
the model (Fig. 4.2). It is this point in the model which is 
used for the thermal control of the heater. Three situations are 
considered for each size container: 

(a) Variation of -the heater cut-off temperature (150F, 155F 
and 160F for a given heater cut-on temperature (140F) as shown in 
Fig, 5.4. 

(b) Variation of the heater cut-on temperature (135F, 140F, 
145F and 150F) for a given heater cut-off temperature (160F) as 
shown in Fig. 5.5. 

(c) Variation in the level of temperature control for a 

specified temperature difference (lOF) as shown in 

Fig. 5.6. 

In comparing the 6 figures, the most significant effect on 
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the heating times occurs when a constant temperature difference 

T - T is maintained, permitting the range of the thermal 
Off on 
control to change. 

5.4 Effect of Initial Temperature 

The initial temperature of the nutrient material is the most 
influential factor on heating-time requirements. Three initial 
temperatures of food substances -are involved, namely 130F, 60F, 
and -lOF. The thermal response of each container is shown in 
Fig. 5.7a and 5.7b for the 60F and 130F initial conditions. The 
thermal response of each container for an initial remperature of 
-lOF is shown in Fig. 5,7c. It should be noted that the predic- 
tions shown in Fig. 5.7c do not include the energy of melting 
associated with the ordinary "thawing" process. The inclusion of 
the phase change increases the heating-time requirement. 

5.5 Effect of Container Size 

Since the size of the containers has already been established 
only the two container sizes to be used in the heating tray are 
included. Typical response characteristics are shown in Fig. 5.7c. 
Parametric studies could be performed to establish the effect of 
length-to-diameter ratios for the cylinder. It is desirable during 
the heating process to have the maximum heat transfer surface area 
per given volume of food. For a cylindrical container heated on 
the sides and bottom this size parameter, the ratio of heated 
surface area to volume, is . 

^ + arRL, ^1 + 2 

^ TTFP L L R 

For the containers used in this study , the smaller container has 
a 30% larger size parameter, indicating that the heating is approxi- 
mately 30% faster. This rough estimate is verified in Fig. 5.7c. 
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5.6 Discuss ion 


Ths parameters involved in the parametric investigations 
were (1) thermal diffusivity of the nutrient material, (2) power 
rating of the heater, (3) control temperatures, (4) initial systems 
temperature, and (5) container dimensions. 

A ±40% change in thermal diffusivity caused a +4 to -lOF 
change in the average temperature after 1 hour for the smaller 
container. Also a -0.4 to + 0.8 hr delay in achieving the minimum 
required average temperature of 140F resulted. The most interesting 
result (Fig. 5.3) shows that a reduction in heat flux by a factor 
of 8 increases heating time only marginally. 

The effect of independently varying the lower or upper heater 
control temperature (Fig. 5.4 and 5.5) have only a small effect on 
heating time. However, as expected, adjusting the level of the 
control while keeping the difference between the upper and lower 
temperatures constant has considerable effect (Fig. 5.6). 

The effect of initial temperature is qualitatively predic- 
able. The lower the initial temperature, the longer the time 
required to heat the nutrient material. The significantly longer 
time required to heat "frozen" foods could represent a problem. 

The above results can be summarized to show the effect of 
initial temperature, normalized thermal diffusivity, and the size 
p 32 rameter, respectively. Figure 5.8 indicates the time required 
for the average temperature to reach 140F when the various parame.ters 
are varied independently from the standard configuration. 
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APPENDIX A 


ADAPTATION OF OLCER’S SOLUTION TO REQUIRED SOLUTION 

'X-- 

Olcer [3] solves the problem of the unsteady temperature dis- 
tribution in a right circular solid cylinder of finite length with 
its entire surface subjected to boundary conditions of the* second 
kind (heat flux) , The three-dimensional transient solution has 
the form (quantities are defined in nomenclature and following) : 
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where 


Fi,F 2 ,F 3 are the fluxes at the bottom, top, and radial surfaces, 
respectively. *■ 

q" is the internal heat source function. 

The dot notation ( ) indicates the derivative with respect to 

time. , 

The eigenvalues, X, , are given by 
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and the prime ( denotes differentiation with respect to the argument 
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12 ! 


Ill" 

12^ 


T.3(r,cp,z,t) = fn(F3) = Tos(r,t) = 


RF; 

2k 


!r 2 2j 


where axial symmetry implies . 

Fi 7^ fn(r) 

Fg fn(r) 

Fa y fn(z) 


In particular, for the heating of food 

Fa = 0 


r. _ r. _ rqol^eater on 
Fi - ^3 - \ 0 /heater off 


(r , Cp, z) = (r , z) 


and for axial symmetry 1=0. Therefore 

Oo = 0 

at „ 


- 


kL ^ 


0, = 0 


o - p 

kR ^ 


Equation (A.l) now becomes 







(r,z,t) = J T^(r,z)rdrdz 


0 L 
” 2 




^ [S ^ 1f]fo } ^ ^ A(l--fflfo } 


+ ^ Y y 

R^L Z L 

m=0 n=0 




2(1 + 6 ) Jo=(m„.R) 


[ J J Jo(W„r)cos(^ + 2^)T.(r,z)rdrdz 


_L 0 
2 




/nn mrzNrqrt'l, ^ 


It is noted that 


/nn nirzN - ' p / nir 

+ ~rr^ = J 1 “=°^ t 


niTz . nrr . nir: 

— - sin y- sin — 


2L niT . nn 

— cos — sxn — 
nn 2 2 
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,v: 




so that for v = 1 


then 


(X) j = XJo(X) , 


R 

I 

0 


rJo(P^r)dr 


1 

— J xjj, (X) <ax 




0 


ulr 

1 P d 

" TT^ I M i 


m 


0 


dX 


U 


m 




When Equations (A-2) , (a- 3) , and (A-4) are combined 


L 

2 R 


T (r , z, t) = — 

R^L 


J. J T^(r,z)rdrdz . 


L 0 
2 


, . 2at \ rqQ-j 

VkL kR Ao ; 




UcL\ 4 


- S) - i(S - im-} 


. 16rr 
+ 2 " 
R L 


L 

2 R 


CO CO 

r r 

m=0 n-0 


_ , » /rnr nnz\ ‘ .2 

J(,(u r)cos(— + — — exp(-aX^ t) 
m \_2 L / mn 

2(1 + 6 ) JoAu R) 

no m 


[ J J T. (r.z)rdrdz 


L ° 
2 


Rd~i 


(A-4) 


(A-5) 



-where 


^ \x 

mn m \ L / * 


where 


U R ^ 0 is the mth root of 


and 


J. '(M R) =0 . 
^ m 


For the first heating period the initial i 


T , (r , z) = ,T . 
1 1 


so that the first integral is equal to and i 
involving vanishes (by A-3) . 


‘''^literature is uniform 


'' bther integral 



APPENDIX B 


COMPUTER PROGRAM DESCRIPTION FOR ANALYTICAL SOLUTION 
DISCUSSED IN SECTION 3.0 

K- 

(For Infinite Cylinder) 


The computer program haa essentially four segments: 

a) Input and Initialization 

b) The Heating Phase 

c) The Insulated Phase 

d) Output 

These four segments can be characterized as follows; 

a) Input and Initialization 

The size (R) of the cylinder, the thermo-physical properties 
of the food (a,c,‘k), the heater flux (qo)f the initial tempera- 
ture (T^) , the temperature ranges desired "^max^ 

time increment (At) are required data. 

b) The Heating Phase 

The time is set to zero, then incremented until the wall 
(the hottest point) temperature reaches the maximum allowable 
temperature. Temperature distributions at selected times are 
stored as the cylinder heats. The final temperature distribu- 
tion becomes the initial temperature distribution of the insu- 
lated phase. 

c) The Insulated Phase 

The time is set to zero> then incremented until the wall 
temperature drops to the minimum temperature. Temperature 
distributions at selected times are stored as the wall tempera- 
ture drops. The cold point temperature is checked to see if 
the food has been heated to the minimum temperature. The final 
temperature distribution becomes the initial temperature distri- 
bution of the next heating phase. 

d) Output 

Desired temperature, profiles, mean temperatures, cold spot 
temperatures, and hot spot temperatures are printed. 

Figure B.l is a simplified flow diagram of. the computer program. 

The letters to the right correspond to the phase described above. 

63 <. 
























« 

APPENDIX C 

COMPUTER PROGRAM DESCRIPTION FOR NUMERICAL SOLUTION 
DISCUSSED IN SECTION 3.0 

(Without Thawing) 

The program follows closely the path outlined in Appendix 
B. Figure C.l as a simplified flow diagram, of the computer 
program. 





